Analysis of aftershocks in a lithospheric model with seismogenic zone governed by damage rheology

We perform analytical and numerical studies of aftershock sequences following abrupt steps of strain in a rheologically layered model of the lithosphere. The model consists of a weak sedimentary layer, over a seismogenic zone governed by a viscoelastic damage rheology, underlain by a viscoelastic upper mantle. The damage rheology accounts for fundamental irreversible aspects of brittle rock deformation and is constrained by laboratory data of fracture and friction experiments. A 1-D version of the viscoelastic damage rheology leads to an exponential analytical solution for aftershock rates. The corresponding solution for a 3-D volume is expected to be sum of exponentials. The exponential solution depends primarily on a material parameter R given by the ratio of timescale for damage increase to timescale for accumulation of gradual inelastic deformation, and to a lesser extent on the initial damage and a threshold strain state for material degradation. The parameter R is also inversely proportional to the degree of seismic coupling across the fault. Simplifying the governing equations leads to a solution following the modified Omori power-law decay with an analytical exponent p = 1. In addition, the results associated with the general exponential expression can be fitted for various values of R with the modified Omori law. The same holds for the decay rates of aftershocks simulated numerically using the 3-D layered lithospheric model. The results indicate that low R values (e.g. R ≤ 1) corresponding to cold brittle material produce long Omori-type aftershock sequences with high event productivity, while high R values (e.g. R ≥ 5) corresponding to hot viscous material produce short diffuse response with low event productivity. The frequency-size statistics of aftershocks simulated in 3-D cases with low R values follow the Gutenberg–Richter power law relation, while events simulated for high R values are concentrated in a narrow magnitude range. Increasing thickness of the weak sedimentary cover produces results that are similar to those associated with higher R values. Increasing the assumed geothermal gradient reduces the depth extent of the simulated earthquakes. The magnitude of the largest simulated aftershocks is compatible with the Båth law for a range of values of a dynamic damage-weakening parameter. The results provide a physical basis for interpreting the main observed features of aftershock sequences in terms of basic structural and material properties.